f{Asfl-C&- nsTW 


(NASA-CR-195799) ANALYSIS OF 

CRUSTAL STRUCTURE OF VENUS 

UTILIZING RESIDUAL L I Nt-OF-S I GHT 

(LOS) GRAVITY ACCELERATION AND 

SURFACE TOPOGRAPHY DATA. A TRIAL OF 

GLOBAL MODELING OF VENUS GRAVITY 

FIELD USING HARMONIC SPLINE METHOD G3/91 

Final Technical Report (Woods Hole 

Oceanographic Inst.) 9 p 


N94-29737 
Unci as 
0004062 


NASA-CR-195799 


FINAL TECHNICAL REPORT 


Analysis of Crustal Structure of Venus Utilizing Residual Line-of-Sight (LOS) 
Gravity Acceleration and Surface Topography Data 

NAG 2-768 

A TRIAL OF GLOBAL MODELING OF VENUS GRAVITY 
FIELD USING HARMONIC SPLINE METHOD 


J./0& 2 — 

°tf 


Ming Fang and Carl Bowin 

Department of Geology &. Geophysics, Woods Hole Oceanographic Institution 

Woods Hole, Ma. 02543 


1. Introduction 

To construct Venus gravity disturbance field ( also referred to as gravity anomaly 
in geological society ) with the spacecraft-obsever line of sight (LOS) acceleration per- 
turbation data, both global approach and local approach can be used. The global 
approach, e.g. spherical harmonic coefficients (Bills et. al. 1987) and local approach, 
e.g. integral operator method (Barriot & Balmino, 1992), based on geodetic techniques 
are generally not the same, so that they must be used separately for mapping long 
wavelength features and short wavelength features. Harmonic spline, as an interpola- 
tion and extrapolation technique, is intrinsically flexible to both global and local map- 
ping of a potential field. Theoretically, It preserves the information of the potential 
field up to the bound by sampling theorem regardless whether it is global or local 
mapping, and is never bothered with truncation errors. A patch by patch construction 
of Venus gravity field at a constant altitude using harmonic spline method has been 
proceeded by Bowin et. al.(1985). In the present investigation we try for the global 
mapping. 

Since the oblateness of Pioneer Venus orbit is extra large sharply reducing the 
magnitude of signals along trajectories away from the periapses, it is unlikely that glo- 
ble modeling with such unevenly magnitudinous data will produce very accurate short 
wavelength gravity features of the planet. Nevertheless, there are other reasons that 
warrant this trial. First, harmonic spline, the technique itself, needs modifications both 
theoretically and n ume rically in order to deal with a fairly large data set of single com- 
ponent measurements, and this has never been done before. Secondly, the only previ- 
ous global modeling of a potential field with harmonic spline was the downward con- 
tinuation of the earth surface magnetic field to the core mantle boundary (Shure et. al. 
1982, Parker & Shure 1982), in which case the data was constantly distributed with 
respect to the radius of the earth. The present study will tell us how the uneven radius 
distribution of data of a harmonic field affects the modeling. This type of information 
will serve as a upper bound of the radial effect for the future processing of the data 
from Magellan mission where the orbit is supposed to be nearly round. 

The emphasis of this report will be on the improvement of harmonic spline 
methodology. In the following sections we describe the new basis functions used in 
this study, then present a singlar value decomposition(SVD) based modification to 
Parker & Shure’s(1982) (hereafter referred to PS) numerical procedure, and finally 
show some of the preliminary results. 
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2. Basis functions for single component incomplete vector data 

Harmonic spline was originally developed for a set of complete vector measure- 
ments of the gradient of a harmonic field which is regular at infinity (Shure ea. al 1982 
hereafter referred to as SPB). At each observation localtion there are three components 
of measurements, and the subset of basis indexed by the location contains three ele- 
ments each of which corresponds to one component of measurement ( see SPB for 
detail). In space geodetic problems, more often we have only one component of meas- 
urement at each observation location like the LOS acceleration perturbation data. 
Under such circumstance the three element sub-basis still works though it may not be 
theoretically optimal. Bowin et. al. (1985) used such basis set to construct their Venus 
gravity map within each patch of area at several constant altitudes. 

It is easy to prove that the optimal representation of the gradient of a harmonic 
field in the sense of minimizing certain norms can be achieved with a basis set that has 
a single element corresponding to each observation location if only a single component 
of measurement for the gradient is available at each observation location. In fact, the 
general theory by SPB does not require that there must be three components at each 
point. Therefore we can proceed strictly in parallel to SPB. Details referred to SPB 
here we simply brief the outcome. Notations except stated are also referred to SPB. 

Let ft j be a unit direction vector at j th observation location with its Cartesian 
components 

( nj- , nf , nf ) 

By regarding each element g ; , • • ■ N) selected from the Hilbert space H 

being the element of a functional for the gradient that gives value of y,- , which is the 
fl j direction component of the gradient of the harmonic field measured at the j th 
location, we have the norm minimizing element B 0 in H such that 

B 0 = Za, g , (1) 

;=i 

and the gradient of the harmonic field B(r) is modeled as 

N 

B(r) = X a ; G ; (r) (2) 

;= i 

where 

G;(r) = v 2 -rrVTr-jt— >' +t i7(8.<t>) 

l/n L + l dfi j { r j J 

For comparison The T matrices in three diferent cases, SPB, Bowin et. al. (1985), and 
this report are listed here in the same order 

T(/, 

3 

r(/, /,;) = I 

i = 1 
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3 3 

T(I,J) = r(/, i , J , j)n}nj 

i= 1;'=1 

where / , J - \ ,2 N represent observation locations, and i , j = 1 , 2 , 3 the 

Cartesian components of vectors. 

Equally important to the optimal nature of the new basis is the fact that the elimi- 
nation of two basis elements for each observation location may save as much as two 
thirds in computer memery for equal number of basis locations without losing, at 
least in theory , resolution power, this can be seen from the expressions of the last two 
r matrices. This tremendous saving of memery makes it possible to use larger number 
of basis locations or even completely non-depleted basis for modeling ( the concept of 
depleted basis is referred to PS ). Roughly speaking, the minimum resolvable 
wavelength by harmonic spline for evenly distributed basis points over a sphere is 
twice the length between each adjacent basis locations( sampling theorem). This 
means that we can increase the resolution power up to a factor of V3* with the same 
computer capacity. 


3. SYD based numerical algorithm 


From numerical point of view, the new basis introduced here is more susceptible 
to rounding errors than the other two. Particularly when directions of the single com- 
ponent data are biased. This is because of that the element corresponding to each 
observation location in the new basis is direction oriented, while for SPB and Bowin 
et al. (1985) bases three elements corresponding to each observation location make 
the bases totally neutral in terms of orientation. Our numerical experiments with 
Pioneer Venus LOS acceleration data using PS’s QR decomposition based algorithm 
has shown that for data distributions of one point at every 6°x6° grid and denser the 
numerical process failed at Cholesky factorization with the non-depleted and very close 
to non-depleted new bases, clearly due to the rounding errors in QR decomposition. 
To solve this problem we modified PS’s algorithm by means of singular value decom- 
position (SVD). 

Without losing generality, let us consider the minimization problem pos^d by PS. 


min 

P 


p r H(5 + YdKP-yl 2 -S 2 ) 


(3) 


IK(3 - yl 2 = S 2 

where X is the Lagrange multiplier, S 2 is the given squared misfit, 

P = ( ft, h f 

is the parameter column to be determined (equivalent to a parameters in eqn. (1) and 


( 2 ) ), 


Y= (Yi.Tfc " • . Yn ) T 


is the observation data column. 


Submitted to NASA meeting 


Hij i=l, 2, • • • , L, 7=1, 2, • • • , L 

are the elements of a positive definite basis matrix H ( depleted or non-depleted), 

Kij i=l, 2 , • • • , V, 7 = 1 , 2 , • • • , L 
are the elements of observation matrix K with the property 

H (K r K) = (K r K) H 

Introducing SVD decomposition of K (e.g. Lawson & Hanson 1974) 

K = U I V r (4) 

where U is a NxN orthogonal matrix, V is a LxL orthogonal matrix, L is a iVxL gen- 
eralized diagonal matrix with none negative diagonal elements ©j, © 2 , ■ ■ • , co L . Fol- 
lowing the logic of PS we can complete from (3) and (4) an iteration process for solv- 
ing P and X. 

(D + i-L 2 )X = ^- (5) 

S 2 = 5 2 in + X 2 X r D 2 X (6) 


^7+1 - ^7 “ 


S 2 Q^j)Sd£sired. . 


35 2 


- y= 0 ,U, • • 


dXj 


Ml 

dX 


= 2X r D 


D+ x £2 


-1 


DX 


Xo = 


C* ■ c 

^ ^ mm 


1 

2 


,d r Z -1 D2“ 2 Dr- 7 'd 
where subscript 7 indicates the iteration times, and 

Smin =/y- d r d 


(7) 

( 8 ) 

(9) 

( 10 ) 


d = U r y 


(ID 


D = Y r HV 


( 12 ) 


X=V r p (13) 

The starting trial value of the Lagrange multiplier X can be calculated from (9) after 
the initial SVD decomposition, then Cholesky factorization or some equivalent method 
should be used to solve for X from eqn. (5) and to perform the iteration until the cal- 
culated misfit falls into the tolerence bound for the desired misfit Finally, (13) 

is solved for p. As a test we repeated, using SVD, the calculations for several patch 
of areas over Venus, previously done with the QR algorithm. For well conditioned 
matrices K, SVD’s results are in good agreement with QR’s as they should. For ill 
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conditioned matrices K, SVD succeeded in one out of total five cases that failed with 
QR. It seems to be a little advantage of SVD over QR, but not much. 

The real advantage appears when the basis is not depleted. In this situation, 

L = N 


K = H 


U = V 


D = I 


where L hence D is strictly diagonal. Diagonaliztion of D greatly simplifies the whole 
process and makes it unneccessary to use the Cholesky factorization or anything simi- 
lar. The new algorithm for a non-depleted basis is as simple as follows. 


di 

X + CDi 


i=l, 2, 


■,N 


(5a) 


S 2 = X 2 X r X 


(6a) 


ax” 


X+co! 


asi 

ax 


= 2XX t 


0 


0 


X+CO, 


w 


Xo = 


cl 

° desired 

w d r Zr 1 d . 


(7a) 


(8a) 


(9a) 


where x t and di denote the i th components of X and d respectively. 

After initial singular value decomposition the rest is just a simple algbera. The 
algorithm of (5a) through (9a) can be analogous to the damping least square algorithm 
with X as the damping factor. Since X has always to be positive (SPB), it will keep 
rounding errors from being amplified by some very small co i s and guarantee the con- 
vergence of the iteration process. In another counterpart test to the one mentioned 
above, four out of five QR failed calculations have been successfully through at very 
fast convergence rate. Aother QR failed case is with a non-depleted basis not applica- 
ble to the new algorithm. In another numerical experiment we picked up the starting 
value of X, Xq , randomly along the positive side of the number axis, and found that 
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iterations nver failed to converge no matter how ill conditioned the matrix K would 
be. 

One should be reminded, however, that a guaranteed convergence by no means 
leads to a guaranteed resolution in the solution. Ill conditioning, as it does to the least 
square solutions, reduces the resolution power by smoothing off short wavelength sig- 
nals. This can be clearly explained by the expression of (5a) and (8a). With extremely 
small co, the denominators are effectively a constant X, and information contained in 
these small co,- is smeared off from the solution. Detailed discussions over ill condi- 
tioning is essentially the same as that for least square problems which can be found in 
from e.g. Lawson & Hanson (1974) and other monographs on least square adjustment. 

Calculations in the rest of this paper are all performed on non-depleted basis with 
the greatly simplified SVD algorithm. 

4. Some preliminary result of global modeling of Venus gravity 

Pioneer Venus orbiter covers Venus surface effectively between -45 and +75 lati- 
tude. Below -45 latitude also has some data coverage but not as regular as the major 
covered region between -45 and 75 latitude. Our so called global modeling is literally 
referred to this major covered region, though the calculation has extended to -80 to 
+80 latitude. Detailed description of the acquisition of the data set has been given in a 
number of previous studies ( Sjogren et.al, 1980, Williams eLal, 1983, Mottinger etal, 
1985, Bills et.al, 1987 and so on ). Reported here are Venus radial gravity anomaly 
and geoid anomaly at 300 km altitude (Fig. 5 and Fig. 6) modeled from LOS accelera- 
tion perturbation data at a reduced density of angular distribution of each grid point 
representing a 4°x4° area. 

Totally we presently are using 2923 points angularly regularly distributed over the 
major covered region. This is a great sacrifice of the data set containing as many as 
145,000 data points, but is the limit that we can handle for the non-depleted basis with 
our present computer facility. The data points are selected from orbits ranging from 
number 111-117, and 418-602 previously used by Bowin eLal (1985) in their patch by 
patch modeling. The radial distribution of data is illustrated by the altitude versus lati- 
tude plot for the orbit number 565 (Fig. 1), and the radial variation from orbit to orbit 
is shown in Fig. 2 where data points of many orbits are plotted. These low altitude 
orbits (150-350 km at periapsis) suffered more from the atmospherical drags than the 
high altitude orbits (950-1350 km at periapsis) used by Mottinger et.al (1985) to form 
their tenth-degree and tenth-order model. On the other hand, low altitude data contain 
more detailed gravity features than the high altitude data. Bills et.al (1987) used low 
altitude data in addition to high altitude one to produce their eighteenth degree and 
eighteenth order model. 

Under perfect condition the shortest wavelength resolvable by the data distribu- 
tion would be near 1000 km equivalent to about 40 spherical harmonic degree. Practi- 
cally this goal can not be achieved due to observation errors and other limitations. 
Aside from the errors in the data itself, there are two effects that will diminish the 
resolution and introduce spurious information. One is the orientation of the LOS data. 
Since the orbital radius of Pioneer Venus is extremely small compare to the distance 
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between Venus and the Earth, the LOS data are very oriented latitudinally. As dis- 
cussed earlier, this worstens the condition of the observation matrix K ( or H identical 
to K ). Aother effect is the radial reduction of the gravitational signal at high altitude. 
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